#dat_full <- read.csv(here("data/dat_04_04_2023.csv"))#dat_full <- read.csv(here("data/dat_28_06_2023.csv"))dat_full <-read.csv(here("data/dat_19_07_2023.csv"))# add phylogenetic tree - only topologies# TODO? - we could get better tree from birdtree.org# we can do 50 different trees as in # https://academic.oup.com/sysbio/article/68/4/632/5267840tree_top <-read.tree(here("R/birds_MA.tre"))# tree with branch lengthstree <-compute.brlen(tree_top)#plot(tree)# turning it into a correlation matrixcor_tree <-vcv(tree,corr=T)
3 Custom functions
Code
# custom functions#' Title: Contrast name generator#'#' @param name: a vector of character stringscont_gen <-function(name) { combination <-combn(name, 2) name_dat <-t(combination) names <-paste(name_dat[, 1], name_dat[, 2], sep ="-")return(names)}#' @title get_pred1: intercept-less model#' @description Function to get CIs (confidence intervals) and PIs (prediction intervals) from rma objects (metafor)#' @param model: rma.mv object #' @param mod: the name of a moderator get_pred1 <-function (model, mod =" ") { name <-firstup(as.character(stringr::str_replace(row.names(model$beta), mod, ""))) len <-length(name)if (len !=1) { newdata <-diag(len) pred <- metafor::predict.rma(model, newmods = newdata,tau2.levels =1:len) }else { pred <- metafor::predict.rma(model) } estimate <- pred$pred lowerCL <- pred$ci.lb upperCL <- pred$ci.ub lowerPR <- pred$cr.lb upperPR <- pred$cr.ub table <-tibble(name =factor(name, levels = name, labels = name), estimate = estimate,lowerCL = lowerCL, upperCL = upperCL,pval = model$pval,lowerPR = lowerPR, upperPR = upperPR)}#' @title get_pred2: normal model#' @description Function to get CIs (confidence intervals) and PIs (prediction intervals) from rma objects (metafor)#' @param model: rma.mv object #' @param mod: the name of a moderator get_pred2 <-function (model, mod =" ") { name <-as.factor(str_replace(row.names(model$beta), paste0("relevel", "\\(", mod,", ref = name","\\)"),"")) len <-length(name)if(len !=1){ newdata <-diag(len) pred <-predict.rma(model, intercept =FALSE, newmods = newdata[,-1]) }else { pred <-predict.rma(model) } estimate <- pred$pred lowerCL <- pred$ci.lb upperCL <- pred$ci.ub lowerPR <- pred$cr.lb upperPR <- pred$cr.ub table <-tibble(name =factor(name, levels = name, labels = name), estimate = estimate,lowerCL = lowerCL, upperCL = upperCL,pval = model$pval,lowerPR = lowerPR, upperPR = upperPR)}#' @title mr_results#' @description Function to put results of meta-regression and its contrasts#' @param res1: data frame 1#' @param res1: data frame 2mr_results <-function(res1, res2) { restuls <-tibble(`Fixed effect`=c(as.character(res1$name), cont_gen(res1$name)),Estimate =c(res1$estimate, res2$estimate),`Lower CI [0.025]`=c(res1$lowerCL, res2$lowerCL),`Upper CI [0.975]`=c(res1$upperCL, res2$upperCL),`P value`=c(res1$pval, res2$pval),`Lower PI [0.025]`=c(res1$lowerPR, res2$lowerPR),`Upper PI [0.975]`=c(res1$upperPR, res2$upperPR), )}#' @title all_models#' @description Function to take all possible models and get their results#' @param model: intercept-less model#' @param mod: the name of a moderator all_models <-function(model, mod =" ", type ="homo") {# getting the level names out level_names <-levels(factor(model$data[[mod]])) dat2 <- model$data mod <- mod VCV1 <-vcalc(vi = dat2$Vd,cluster = dat2$SubjectID,rho =0.5)#model$data[[mod]] <- factor(model$data[[mod]], ordered = FALSE)# meta-regression: contrasts # helper function to run metafor meta-regression run_rma1 <-function(name) {rma.mv(yi = SMD, V = VCV1, mods =~relevel(dat2[[mod]], ref = name), random =list(~1| FocalSpL,~1| RecNo,~1| Obs_ID),test ="t",method ="REML",sparse =TRUE,data = dat2) } run_rma2 <-function(name) {rma.mv(yi = SMD,V = VCV1,mods =~relevel(dat2[[mod]], ref = name),random =list(~1| FocalSpL,~1| RecNo,formula(paste("~", mod, "| Obs_ID"))),test ="t",method ="REML",sparse =TRUE,data = dat2) }# results of meta-regression including all contrast results; taking the last level out ([-length(level_names)])# this does not work for hetero model?if (type =="homo"){ model_all <-map(level_names[-length(level_names)], run_rma1) } else { model_all <-map(level_names[-length(level_names)], run_rma2) }# getting estimates from intercept-less models (means for all the groups) res1 <-get_pred1(model, mod = mod)# getting estiamtes from all contrast models res2_pre <-map(model_all, ~get_pred2(.x, mod = mod))# a list of the numbers to take out unnecessary contrasts contra_list <-Map(seq, from=1, to=1:(length(level_names) -1)) res2 <-map2_dfr(res2_pre, contra_list, ~.x[-(.y), ]) # creating a table res_tab <-mr_results(res1, res2) %>%kable("html", digits =3) %>%kable_styling("striped", position ="left") %>%scroll_box(width ="100%")# results res_tab}# test #all_models(mod1, mod = "Treat_mod")#all_models(mod1s, mod = "Treat_mod")#all_models(mod1s2, mod = "Treat_mod")
4 Preparing data set
Code
# function for calculating varianceVd_func <-function(d, n1, n2, design, r =0.5){# independent designif(design =="among"){ var <- (n1 + n2) / (n1*n2) + d^2/ (2* (n1 + n2 -2)) # variance } else { # dependent design var <-2*(1-r) / n1 + d^2/ (2*(n1 -1)) # variance } var # return variance}# getting Hedges' g - get small size bias corrected effect sizedat_full$SMD <- dat_full$d / (1-3/(4* (dat_full$NTreat + dat_full$Ncontrol) -9))# flipping d dat_full$SMD <- dat_full$d*dat_full$Direction*dat_full$PredictedDirection# calucating Vddat_full$Vd <-with(dat_full, pmap_dbl(list(SMD, NTreat, Ncontrol, Design), Vd_func))# extra useful function# function for getting mean and sd from median, quartiles and sample size# get_mean_sd <- function(median, q1, q3, n){# sd <- (q3 - q1) / (2 * (qnorm((0.75 * n - 0.125) / (n + 0.25)))) # sd# mean <- (median + q1 + q3)/3 # mean# c(mean, sd)# }# observation iddat_full$Obs_ID <-1:nrow(dat_full)dat_full$Phylo <-gsub(" ", "_", dat_full$FocalSpL)# filtering very large variance and also very small sample sizedat_int <- dat_full %>%filter(Vd <10& Ncontrol >2& NTreat >2)#dim(dat_full)#dim(dat_int)# sorting out modality stuff# creat - 1,2,3 modality - also easier classification A, O, V (AOV = L) dat_int %>%mutate(Treat_mod =case_when(Treatment =="A"~"A", Treatment =="AV"~"AV", Treatment =="AVG"~"AV", Treatment =="AVM"~"AV", Treatment =="L"~"AVO", Treatment =="O"~"O", Treatment =="OV"~"OV", Treatment =="V"~"V", Treatment =="VG"~"V", Treatment =="VM"~"V", Treatment =="VP"~"V"),# into how manyTreat_No =case_when(Treatment =="A"~1, Treatment =="AV"~2, Treatment =="AVG"~2, Treatment =="AVM"~2, Treatment =="L"~3, Treatment =="O"~1, Treatment =="OV"~2, Treatment =="V"~1, Treatment =="VG"~1, Treatment =="VM"~1, Treatment =="VP"~1),# des it have some add-onsAdd_on =case_when(Treatment =="A"~"No", Treatment =="AV"~"No", Treatment =="AVG"~"Yes", Treatment =="AVM"~"Yes", Treatment =="L"~"No", Treatment =="O"~"No", Treatment =="OV"~"No", Treatment =="V"~"No", Treatment =="VG"~"Yes", Treatment =="VM"~"Yes", Treatment =="VP"~"Yes"), ) -> dat# creating data just for A, V, and AV dat_short <- dat %>%filter(Treat_mod =="A"| Treat_mod =="V"| Treat_mod =="AV")# for add-on, we only need V and AVdat_short_add <- dat %>%filter(Treat_mod =="AV"| Treat_mod =="V")dat <- dat %>%mutate_if(is.character, as.factor)# reordering factors for better visualization# dat$Treat_mod <- factor(dat$Treat_mod, levels = c("A", "V", "AV", "O", "OV", "AVO"))
5 Exploratory visualization
For Treat_mod (Treatment), we will only visualise A, V, and AV as O $r $, OV, and AVO are much rarer. But for Type (Trait type), we will use all data.
orchard_plot(mod_trt_2, mod ="Treat_mod",group ="RecNo", xlab ="Standardised mean differnece (SMD)",branch.size =3)
Code
# result table (get all the constrasts)all_models(mod_trt_2, mod ="Treat_mod")
Fixed effect
Estimate
Lower CI [0.025]
Upper CI [0.975]
P value
Lower PI [0.025]
Upper PI [0.975]
A
0.410
0.211
0.609
0.000
-1.497
2.317
V
0.440
0.223
0.658
0.000
-1.578
2.458
AV
0.448
0.207
0.690
0.000
-1.003
1.900
A-V
0.033
-0.244
0.310
0.815
-1.843
1.909
A-AV
0.073
-0.254
0.399
0.662
-1.811
1.957
V-AV
0.040
-0.267
0.346
0.800
-1.841
1.920
Code
# comparison modelsanova(mod_trt_1, mod_trt_2) # modeling heteroscedasticity is better
df AIC BIC AICc logLik LRT pval QE
Full 8 1692.5188 1727.6541 1692.7637 -838.2594 NA
Reduced 6 1707.4108 1733.7623 1707.5532 -847.7054 18.8921 <.0001 NA
Code
# the effect of additions# this is a part of sensitivity analysisdat_short_add$Treat_add <-as.factor(paste(dat_short_add$Treat_mod , dat_short_add$Add_on, sep ="-"))VCVs2 <-vcalc(vi = dat_short_add$Vd,cluster = dat_short_add$SubjectID,rho =0.5)mod_add <-rma.mv(yi = SMD, V = VCVs2, random =list(~1|FocalSpL , ~1| RecNo, ~1| Obs_ID), mod =~ Treat_add -1,test ="t",struct ="DIAG",method ="REML", sparse =TRUE,data = dat_short_add)summary(mod_add)
orchard_plot(mod_add, mod ="Treat_add",group ="RecNo", xlab ="Standardised mean differnece (SMD)",branch.size =3)
Code
# result table (get all the constrasts)all_models(mod_add, mod ="Treat_add")
Fixed effect
Estimate
Lower CI [0.025]
Upper CI [0.975]
P value
Lower PI [0.025]
Upper PI [0.975]
AV-No
0.427
0.150
0.704
0.003
-1.366
2.220
AV-Yes
0.572
0.015
1.128
0.044
-1.285
2.428
V-No
0.426
0.211
0.642
0.000
-1.358
2.210
V-Yes
0.293
-0.118
0.704
0.162
-1.525
2.111
AV-No-AV-Yes
0.145
-0.466
0.755
0.642
-1.729
2.018
AV-No-V-No
-0.001
-0.324
0.322
0.996
-1.801
1.800
AV-No-V-Yes
-0.134
-0.619
0.351
0.587
-1.970
1.702
AV-Yes-V-No
-0.145
-0.734
0.443
0.627
-2.012
1.721
AV-Yes-V-Yes
-0.279
-0.936
0.379
0.405
-2.168
1.611
V-No-V-Yes
-0.133
-0.578
0.312
0.556
-1.959
1.693
Code
# testing the number of stimulimod_ord <-rma.mv(yi = SMD, V = VCV, random =list(~1|FocalSpL , ~1| RecNo, ~1| Obs_ID), mod =~ Treat_No, test ="t",method ="REML", sparse =TRUE,data = dat)summary(mod_ord)
bubble_plot(mod_ord,mod ="Treat_No",group ="RecNo",xlab ="The number of simuli",g =TRUE)
Code
#homoscedastic model# Type of responses mod_type1 <-rma.mv(yi = SMD, V = VCV, random =list(~1| FocalSpL,~1| RecNo,~1| Obs_ID),mod =~ Type -1,test ="t",method ="REML",sparse =TRUE,data = dat)summary(mod_type1)
orchard_plot(mod_type2, mod ="Type",group ="RecNo",xlab ="Standardised mean differnece (SMD)",branch.size =3)
Code
# result tableall_models(mod_type2, mod ="Type")
Fixed effect
Estimate
Lower CI [0.025]
Upper CI [0.975]
P value
Lower PI [0.025]
Upper PI [0.975]
Behaviour
0.503
0.361
0.644
0.000
-1.275
2.281
LifeHistory
0.353
0.147
0.559
0.001
-1.015
1.721
Physiology
0.012
-0.292
0.316
0.938
-2.328
2.352
Behaviour-LifeHistory
-0.178
-0.418
0.062
0.146
-1.975
1.619
Behaviour-Physiology
-0.458
-0.723
-0.193
0.001
-2.259
1.342
LifeHistory-Physiology
-0.280
-0.586
0.025
0.072
-2.087
1.526
Code
# heteroscadasticity model better than the homoscedasticity model# note LifeHistory has lowest variation but this may be expected? # as it is less flexiable (e.g. the number of eggs?)anova(mod_type1, mod_type2)
df AIC BIC AICc logLik LRT pval QE
Full 8 1757.5731 1793.2272 1757.8024 -870.7865 NA
Reduced 6 1785.4590 1812.1996 1785.5923 -886.7295 31.8859 <.0001 NA
res<-mod_results(mod_dur_int, mod ="ln_duration", group ="RecNo", by ="Type")bubble_plot(res,mod ="ln_duration",group ="RecNo",xlab ="log(duration in days)",condition.nrow =3,g =TRUE)
Code
# sexmod_sex1 <-rma.mv(yi = SMD, V = VCV, random =list(~1|FocalSpL , ~1| RecNo,~1| Obs_ID),mods =~ Sex -1,test ="t",method ="REML",sparse =TRUE,data = dat)summary(mod_sex1)
orchard_plot(mod_full,mod ="Type",group ="RecNo", xlab ="Standardised mean differnece (SMD)",branch.size =3)
Code
orchard_plot(mod_full,mod ="Treat_mod",group ="RecNo", xlab ="Standardised mean differnece (SMD)",branch.size =3)
Code
# interaction plot of duration with trait typesint_type <-mod_results(mod_full, mod ="sln_duration", group ="RecNo", weights ="prop",by ="Type")bubble_plot(int_type, group ="RecNo", mod ="sln_duration", xlab ="ln(duration in days) (z-transformed)",legend.pos ="top.left", condition.nrow =3, g =TRUE)
Code
# interaction plot of duration with trait typesint_trt <-mod_results(mod_full, mod ="sln_duration", group ="RecNo", weights ="prop",by ="Treat_mod")bubble_plot(int_trt, group ="RecNo", mod ="sln_duration", xlab ="ln(duration in days)",legend.pos ="top.left", condition.nrow =3, g =TRUE)
# displays delta AICc <2candidates_aic2 <-subset(candidates, delta <5) # model averagingmr_averaged_aic2 <-summary(model.avg(candidates, delta <5)) # relative importance of each predictor for all the modelsimportance <-sw(candidates)